Accessing Cloud Optimized GeoTIFF (COG) - S3 Direct Access

Summary

Harmonized Landsat Sentinel-2 (HLS) Operational Land Imager Surface Reflectance and TOA Brightness Daily Global 30m v2.0 (L30) (10.5067/HLS/HLSL30.002)

Requirements

AWS intance running in us-west 2

Earthdata Login

.netrc file

Learning Objectives

S3 Direct Access

Cloud Optimized GeoTIFF (COG)

Using Harmonized Landsat Sentinel-2 (HLS) version 2.0

Import Packages

import os
import requests 
import boto3
from osgeo import gdal
import rasterio as rio
from rasterio.session import AWSSession
import rioxarray
import hvplot.xarray
import holoviews as hv

Single File S3 Direct Access

Temporary Credentials

s3_cred_endpoint = {
    'podaac':'https://archive.podaac.earthdata.nasa.gov/s3credentials',
    'gesdisc': 'https://data.gesdisc.earthdata.nasa.gov/s3credentials',
    'lpdaac':'https://data.lpdaac.earthdatacloud.nasa.gov/s3credentials',
    'ornldaac': 'https://data.ornldaac.earthdata.nasa.gov/s3credentials',
    'ghrcdaac': 'https://data.ghrc.earthdata.nasa.gov/s3credentials'
}

def get_temp_creds(provider):
    return requests.get(s3_cred_endpoint[provider]).json()

temp_creds_req = get_temp_creds('lpdaac')
#temp_creds_req

session = boto3.Session(aws_access_key_id=temp_creds_req['accessKeyId'], 
                        aws_secret_access_key=temp_creds_req['secretAccessKey'],
                        aws_session_token=temp_creds_req['sessionToken'],
                        region_name='us-west-2')

GDAL Configurations

GDAL is a foundational piece of geospatial software that is leveraged by several popular open-source, and closed, geospatial software. The rasterio package is no exception. Rasterio leverages GDAL to, among other things, read and write raster data files, e.g., GeoTIFFs/Cloud Optimized GeoTIFFs. To read remote files, i.e., files/objects stored in the cloud, GDAL uses its Virtual File System API. In a perfect world, one would be able to point a Virtual File System (there are several) at a remote data asset and have the asset retrieved, but that is not always the case. GDAL has a host of configurations/environmental variables that adjust its behavior to, for example, make a request more performant or to pass AWS credentials to the distribution system. Below, we’ll identify the evironmental variables that will help us get our data from cloud

rio_env = rio.Env(AWSSession(session),
                  GDAL_DISABLE_READDIR_ON_OPEN='TRUE',
                  GDAL_HTTP_COOKIEFILE=os.path.expanduser('~/cookies.txt'),
                  GDAL_HTTP_COOKIEJAR=os.path.expanduser('~/cookies.txt'))
rio_env.__enter__()

<rasterio.env.Env at 0x7f5e4b7fa490>

s3_url = 's3://lp-prod-protected/HLSL30.020/HLS.L30.T11SQA.2021333T181532.v2.0/HLS.L30.T11SQA.2021333T181532.v2.0.B04.tif'

da = rioxarray.open_rasterio(s3_url)

da

da_red = da.squeeze('band', drop=True)
da_red

<xarray.DataArray (y: 3660, x: 3660)>
[13395600 values with dtype=int16]
Coordinates:
  * x            (x) float64 7e+05 7e+05 7e+05 ... 8.097e+05 8.097e+05 8.097e+05
  * y            (y) float64 4.1e+06 4.1e+06 4.1e+06 ... 3.99e+06 3.99e+06
    spatial_ref  int64 0
Attributes:
    _FillValue:    -9999.0
    scale_factor:  0.0001
    add_offset:    0.0
    long_name:     Red

xarray.DataArray

y: 3660
x: 3660

...
```
[13395600 values with dtype=int16]
```
Coordinates: (3)
- x
  (x)
  float64
  7e+05 7e+05 ... 8.097e+05 8.097e+05
```
array([699975., 700005., 700035., ..., 809685., 809715., 809745.])
```
- y
  (y)
  float64
  4.1e+06 4.1e+06 ... 3.99e+06
```
array([4100025., 4099995., 4099965., ..., 3990315., 3990285., 3990255.])
```
- spatial_ref
  ()
  int64
  0
  crs_wkt :
  PROJCS["UTM Zone 11, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
  semi_major_axis :
  6378137.0
  semi_minor_axis :
  6356752.314245179
  inverse_flattening :
  298.257223563
  reference_ellipsoid_name :
  WGS 84
  longitude_of_prime_meridian :
  0.0
  prime_meridian_name :
  Greenwich
  geographic_crs_name :
  Unknown datum based upon the WGS 84 ellipsoid
  horizontal_datum_name :
  Not_specified_based_on_WGS_84_spheroid
  projected_crs_name :
  UTM Zone 11, Northern Hemisphere
  grid_mapping_name :
  transverse_mercator
  latitude_of_projection_origin :
  0.0
  longitude_of_central_meridian :
  -117.0
  false_easting :
  500000.0
  false_northing :
  0.0
  scale_factor_at_central_meridian :
  0.9996
  spatial_ref :
  PROJCS["UTM Zone 11, Northern Hemisphere",GEOGCS["Unknown datum based upon the WGS 84 ellipsoid",DATUM["Not_specified_based_on_WGS_84_spheroid",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-117],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["Easting",EAST],AXIS["Northing",NORTH]]
  GeoTransform :
  699960.0 30.0 0.0 4100040.0 0.0 -30.0
```
array(0)
```
Attributes: (4)
_FillValue :
-9999.0
scale_factor :
0.0001
add_offset :
0.0
long_name :
Red

da_red.hvplot.image(x='x', y='y', cmap='gray', aspect='equal')

Unable to display output for mime type(s):

rio_env.__exit__()